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Abstract. - The viscosity and self-diffusion constant of a mesoscale hydrodynamic method, 
dissipative particle dynamics (DPD), are investigated. The viscosity of DPD with finite time 
step, including the Lowe- Anderson thermostat, is derived analytically for the ideal-gas equation 
of state and phenomenologically for systems with soft repulsive potentials. The results agree well 
with numerical data. The scaling of the local relative velocity in molecular dynamics simulations 
is shown to be useful to obtain faster diffusion than for the DPD thermostat. 



Introduction. — Soft matter systems such as poly- 
mer solutions, colloidal suspensions, vesicles, cells, and 
microemulsions exhibit many interesting dynamical be- 
haviors, where hydrodynamic flow plays an important 
role, as do thermal fluctuations. The characteristic length 
(ran to fim) and time (ns to s) scales of soft-matter sys- 
tems are typically much larger than the atomistic scales. 
Coarse-grained molecular models and simulation methods 
are therefore necessary to simulate mesoscale phenomena 
with reasonable computational effort. Dissipative particle 
dynamics (DPD) [1-14] has been developed for this pur- 
pose, and has been applied to various systems such as col- 
loids [6] polymers [2,7,8] and lipid membranes [9], DPD is 
an off-lattice hydrodynamic method, which has two main 
features: soft-repulsive interaction potentials and the pair- 
wise version of a Langevin thermostat. Since there is no 
impenetrable exclude volumes, a DPD particle describes 
not a solvent molecule but a fluid clement, which rep- 
resents clusters of solvent molecules. The main motiva- 
tion for the use of soft potentials is that they allows large 
time steps for the time evolution; however, it has been 
shown that the simulations have to be checked carefully 
in each case by monitoring the configurational tempera- 
ture [15], in order to avoid artifacts due to too large time 
steps [5] . DPD shares many properties with direct simula- 
tion Monte Carlo [16] and multi-particle collision (MPC) 
dynamics [17-21] as pointed out very recently in Ref. [22]. 

The transport coefficients of DPD have been studied 
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for about a decade [1, 10-14]. However, the viscosity of 
DPD has been derived analytically only for the ideal-gas 
equation of state in the small time-step limit [10], and par- 
tially for the Anderson-thermostat [23] version proposed 
by Lowe [24] . Recently, numerical integrators for the DPD 
thermostat [3,4], which have no time-step dependence on 
thermodynamic properties, were proposed; they include 
the Lowe- Anderson thermostat (Lowe-AT) as a specific 
limit. However, the transport coefficients do depend on 
the time step. Therefore, a detailed understanding of the 
transport coefficients for finite time step is very important 
to control and tune the hydrodynamic properties of DPD 
fluids. 

In this letter, we calculate the time-step dependence of 
the viscosity, and in particular investigate the contribution 
due to the interaction potential. The viscosity of DPD 
consists of three contributions, r\ = r\Y m + f]oa\ + %>ot- The 
kinetic viscosity 77km, collision viscosity ?7 co i, and poten- 
tial viscosity 77 po t result from the momentum transfer due 
to particle displacements, collisions generated the DPD 
thermostat (arising from frictional interactions and ther- 
mal noise), and potential interactions, respectively. We 
determine these three contributions both analytically and 
numerically. In previous DPD simulations with a repul- 
sive potential, the contributions of the potential interac- 
tions were often neglected in the discussion of transport 
coefficients. However, the potential contributes to the vis- 
cosity as well as the DPD thermostat in typical simula- 
tion conditions. We also study the self-diffusion constant 
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D of a DPD particle and the ratio of momentum to mass 
transport, which is characterized by the Schmidt number 
Sc = v/D, where v = r)/p is the kinematic viscosity. Fi- 
nally we show that faster relaxation and larger diffusion 
constants D can be obtained in Molecular Dynamics (MD) 
simulations by the rescaling of the local relative velocity 
to control temperature instead of a DPD thermostat. 

Methods. — The DPD thermostat is a modified 
Langevin thermostat, where the friction and noise terms 
are applied to the relative velocities of the neighbor pairs. 
The equation of motion for the i-th particle with mass m 
is given by 




where = Vj - Vj, = r l - rj, r io = [r^-], and 
Tij = Tij/rij. The Gaussian white noise £ij(t) obeys the 
fluctuation-dissipation theorem, with average (i)) = 
and variance (£ij(t)£i'j'(tf)) = 2k B T(5 li >5 jj > + S lj ^S i: j>)S(t- 
t'), where k B T is the thermal energy. This thermostat 
is applied only in the direction to conserve the local 
angular momentum. In DPD, a linear weight function 
w(rij) = wi(r i0 ), with 

U,l(rij)=7 ( 1 ~2;)' (2) 

is typically employed, which vanishes beyond the cutoff at 
r ij = r cut- Furthermore, DPD is usually combined with a 
soft repulsive potential [2], 

with the same cutoff r cut , but other potentials are also 
available. 

The DPD equation (1) is discretized by the Shardlow's 
SI splitting algorithm [3], where each thermostat of the ij 
pair is separately integrated, 

<° w = v * + {-A( r ij)vij ■ Uj + B{r l3 % hn }r tj , 
vf W = Wj - {-A(r tf )v« • % + B(r y )6j,n}?«, (4) 

with 

_ w{r l0 )At/m _ y/wjn^At/m 

{ij> l + w( rij )At/m' l + w( rij )At/m- 

(5) 

The discretized Gaussian noise £y )n is determined by 
(£ij,n£i'j',n') = 2kBT(8u>6jj> + 5ij> 5ij>)S nn > • This split- 
ting algorithm belongs to the generalized Lowe- AT [4], 
because the factors Afcj) and Bfaj) satisfy the relation 
B = y/A(l - A)/m [22]. Thus, for U = this algorithm 
yields the flat radial distribution function of an ideal gas 



for any time step At, without any deviation of the ki- 
netic temperature from the thermostat temperature. In 
the Lowe- AT [24] , the relative velocity v,j f y of a neighbor 
pair ij with rij/r cut < 1 is updated by assigning a random 
number drawn from the Maxwcll-Boltzmann distribution 
with the probability V at each time step At (i.e. velocities 
are updated with the rate T = V /At). When a piecewise 
constant weight function w(rij) = wo(^ij), where 

W ^ = { otherwise, (6) 

is employed, Eq. (4) with ^At/m = 1 gives the Lowe- AT 
for T' = 1. 

The viscosities are calculated from simulations of sim- 
ple shear flow in three dimensions with Lees-Edwards 
boundary conditions [25]. We use the weight w\(rij), 
defined in Eq. (2), and the splitting algorithm (4) for 
the DPD simulations. However, the derived analytical 
expressions can be applied to other weights w(rij) and 
other generalized Lowe- AT algorithms such as A(rtj) in 
table I of Rcf. [4]. The self-diffusion constant D is cal- 
culated from the mean square displacement of a particle, 
({rj(t) — r^O)} 2 ) = 2dDt, where d is the spatial dimension. 

We have performed simulations with the usual soft po- 
tential (3) in order to investigate the effect of interaction 
potentials. The multi-time-step algorithm [4, 26] is em- 
ployed, with a shorter time step 5t for the force —dU/dri, 
so that the configurational [5,15] and thermostat temper- 
atures differ by less than 0.5%. The side lengths of the 
simulation box are L y > 40r cut , L x = L z — 10r cut and 
L x = L y = L z = 20r cut for the calculation of the viscosity 
and the diffusion constant, respectively. The error bars 
of the simulation results are estimated from three inde- 
pendent runs. We display our simulation results in form 
of dimcnsionlcss quantities, indicated by a superscript, 
7* = jTo/m, At* — At/ To, 5t* = St/rQ, and the num- 
ber density n* — nr cut d , which corresponds to measuring 
length, time, viscosity, and diffusion constant of a particle 
in units of r cut , r = r cut \Jm/k B T, rj = VrnksT /rcut^ 1 , 
and Do — r cut \Jk B T /m, respectively. 

Viscosity of an ideal DPD gas. — First, we derive 
the viscosity of an ideal gas of DPD particles, with U = 
and ?7pot = 0, using a kinetic-theory approach. In simple 
shear flow with flow velocity v = jye Xl the xy component 
of the stress tensor is given by a xy = 777. The viscosities 
?7kin and ?7 co i are calculated from the stress due to the 
kinetic and collisional contributions, respectively. 

The kinetic stress cr^ n is the momentum flux due to par- 
ticles crossing a plane of constant y; it can be calculated by 
following the derivation for MPC in Ref. [19]. The stress 
is written as 

< = f ^ i dvv x P(y-iye x ) 

At I J Vy> -^ 

+ dy dv v x P(v-iye x )\, (7) 

Jn J„..s- JL- I 
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Fig. 1: (Color online) Dependence of the viscosity r\ of an ideal 
DPD gas (with U = 0) on (a) the time step At* at n* = 3 
and (b) the number density n for At* = 1. Symbols indicate 
simulation data for 7* = 1 (A, o) and 7* = 9 (o, □). Lines 
represent the analytical results of Eqs. (8) and (12). The inset 
in (a) shows the dependence of the diffusion constant D of the 
ideal DPD gas on the time step At* for n* — 3 and 7* = 9. 



where P(v) is the velocity probability distribution in 
the local rest frame. This stress can be rewritten as 
°asy 1 — rnn^At^) /2 — (v x v y )). The velocity distribu- 
tion is shifted by particle streaming in the time interval 
At, so that (v x v y ) — ► (v x v y ) — jAt(v x v y ). Then, the DPD 
collisions of Eq. (4) modify it as (v x v y ) — > s(v x v y ). Thus, 
the self-consistency condition of a stationary shear flow is 
{v x Vy) = s((v x v y ) — A /At(v x v y )). The kinetic viscosity ?7kin 
is then given by [19] 



nk B TAt 



1 



1- s 



(8) 



The remaining task is to calculate the factor s for the DPD 
collisions. The i-th particle collides with a multitude of 
other particles at the same time step, so that s = (HjSij). 
Eq. (4) together with a molecular chaos assumption im- 
plies Sij = 1 - A(x$j + yfj) + AA 2 xf J yf J , where Xij and 
jjij are the components of Pjj. In an ideal gas, the local 
number density fluctuates around the average n, and the 
number of particles k per volume AV is given by the Pois 




Fig. 2: (Color online) Dependence of the viscosity of an ideal 
DPD gas with Lowe- Anderson thermostat on the normalized 
collision frequency T. Symbols represent simulation data at 
At* = 0.1 (o, □) and At* = 1 (A, o). Lines indicate the 
analytical results. 



the factor s is given by 

s = exp in [ (-2A{r)x 2 + AA{rfx 2 y 2 ) dV \ (9) 



d{d + 2) 



Eqs. (8) and (9) give the kinetic viscosity ?7kin for a finite 
time step At. In the continuum limit At <C 1, we recover 
the result 



dmk^T 
2\w]„ ' 



w. 



g{r)w{r)dV (10) 



of Ref. [10], where g(r) is the radial distribution function, 
with g(r) = 1 for the ideal gas. In the Lowe- AT, the factor 
s is given by s = exp(— 7r&nr') with 6=1/2 and b = 16/45 
in two and three spatial dimensions, respectively. In the 
limit At <C 1 with finite P, 77km = k^T/iibY. 

The collisional stress cr™ 1 is the momentum flux due to 



DPD collisions 
plane at y = yo 



col 2 

o X y = ~n 



determined by Eq. (4) — crossing a 



0. 



dy, 



Vij>Vi 



After substitution of Eq. (4) and 



At 



IVij 



(11) 



into 



Eq. (11) and interchange of the order of integration, 77 co i 
is found to be 



Vcol 



n" l j A(r)mr 2 x 2 y 2 
~ 1 At 



wr 



2d(d+2) 



1 + wAt/m 



(12) 



2 \ in the 



iAy 



son distribution, P(k) = 

(c fc ) = exp{(— 1 + c)nAV} for some constant c. Therefore, 



Equation (12) gives ?7 co i = {n 2 /2d(d + 2)}[wr 
limit At < 1. For the Lowe-AT, Eq. (12) implies 
tcoi = 7rmri 2 rr cut 4 /64 and r] co \ = TTmn 2 Tr cut 5 /75 in two 
and three spatial dimensions, respectively. These results 
(nAV) k /k\, which implies agree with the collisional viscosities obtained in Refs. [10] 

and [24]. 
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The analytical results agree well with the numerical 
data, as shown in Figs. 1 and 2. As At increases, Ty^m in- 
creases but r] co \ decreases, just like the viscosities of MPC- 
Langevin dynamics [22] . Although is almost indepen- 
dent of the density n at small At, 77ki n increases with n 
at large At, see Fig. 1(b). There are small deviations be- 
tween analytical and numerical results in Fig. 1. They 
are of the same order of magnitude as the deviations for 
At <C 1 reported in Ref. [12], which have been explained 
by correlation effects between DPD collisions [12]. In the 
Lowe-AT, the viscosities depend on At for large T, see 
Fig. 2. For V = TAt > 1, our theory overestimates ?7 co i, 
since the relative velocities of some ij pairs are updated 
more than once in one time step. 

Viscosity with interaction potential. With in- 
teraction potential, an additional momentum flux crossing 
a plane at y — 0, is caused by the forces /(ry) = —dU/drij 
between ij pairs with yi > and yj < 0. The potential 
viscosity 77 pot is given by 



»7 P ot 



n 

7 Jn 

2 

n 



dyi 



llr :, 9(nj)f(rij)xij 



Vij >Vi 



-g-r / dV 9(r)f(r)xy, 



(13) 



which is the potential term of the Irving-Kirkwood for- 
mula of the viscosity [27] . The potential also modifies rykin 
with an additional velocity relaxation, while ?7 co i can be 
calculated by Eq. (12) with non-uniform g(r). The viscos- 
ity with an interaction potential has been derived analyti- 
cally for some cases [28], but is generally very complicated. 
Therefore, we employ a simple phenomcnological expres- 
sion instead, and focus on the explanation of qualitative 
dependences. 

In the molecular-chaos approximation, the velocity 
auto-correlation function of a particle in a gas shows an 
exponential decay, (v(t)v(Q)} — exp(— </>£). This behav- 
ior corresponds to the assumption of a Langevin equation; 
dv/dt — —<f)v + ^/(fi^(t)/ra for the particle velocity v in the 
local rest frame. For At <C 1, the DPD collisions generate 
an auto-correlation function with an initial exponential 
decay (for small times t) with </>dpd = n[w] g /dm. Here, 
we assume that the potential also generates an exponen- 
tial auto-correlation function with rate (/> po t , although the 
auto-correlation function determined numerically is not 
exponential, and shows larger deviation from an expo- 
nential decay for larger potential strengths a or particle 
densities n. Then, the kinetic viscosity is given by 



nk^T 



2(0 P ot + 4> 



DPD J 



(14) 



compare Eq. (10). 

In order to estimate ?7pot, an expression for the corre- 
lations of ij pairs is required. We mimic the potential 
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Fig. 3: (Color online) Dependence of the viscosity of DPD with 
soft potentials, Eq. (3), on (a) friction coefficient 7* for n* = 1 
and a — 25, (b) potential strength a for 7* = 4.5, and (c) 
number density n for 7* = 4.5. In all cases, St* — 0.01, and 
At* =0.1. Symbols indicate simulation data for (b) n* — 3 
(o, □, A) and n* = 1 (o, x, +), and (c) a = 5 (o, □, A) and 
a = 25 (o, x, +). Lines for 77 co i represent the analytical results 
of Eq. (12). Lines for r^kin show the results of Eq. (14) with 
4> P ot fitted by Eq. (16). Lines for rj po t are guides to the eye. 



contribution by a DPD thermostat 
rfv 



-jf = { _ 7pot|/(nj)|vij • hj + <r£ij(t)} hj, (15) 



where a = ^7pot 1/(^)1, since the restoring force should 
be proportional to |/(^j)|, and mj in the direction 
Yij. Following the derivation of Eq. (12) with pot = 
n l -pot\\ f \\ g / dm, we obtain the viscosity 



Vpot 



I0pot[|/|? 



(16) 



2(<*+2)[|/|] fl ' 
Thus, the viscosities ?7kin and rj po t can be estimated by 
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Eqs. (14) and (16) with the parameter <^ pot and the radial 
distribution function g(r). 

Figure 3 shows the viscosities of the DPD fluid with 
interaction potential. We calculate g{r) from equilibrium 
simulations, fit (f> po t to rj pot , and then estimate from 
Eq. (14). This underestimates 7/kinj but reproduces very 
well the qualitative dependence on friction coefficient 7, 
potential strength a, and number density n. The kinetic 
viscosity 7/kin decreases with increasing 7 or a. The poten- 
tial viscosity ?7p t is almost independent of 7 and increases 
with a. The collision viscosity ?7coi is almost independent 
of a and shows very good agreement between the theory 
and simulations. 

Diffusion. — Next, we derive the self-diffusion con- 
stant D of an ideal gas of DPD particles (with [7 = 0) for 
finite time steps. Following the derivation of Eq. (9), we 
find that the velocity correlation for one step is given by 
(v x (t + At)v x (t)) = exp(—n[A] g /d). Under the molecular 
chaos assumption, i.e. (v x (kAt)v x (0)) = (v x (At)v x (0)) h , 
the diffusion constant is then given by 

fceTAi / 1 1\ 

m \l-exp(-n[A] g /d) 2 J K ' 

In the limit At <C 1, the diffusion constant becomes 
D = dkBT/n[w] g , in agreement with the result of 
Ref. [10]. However, the velocity auto-correlation func- 
tion (v x (kAt)v x (0)} with large dimensionless friction co- 
efficient 7* has a long-time tail due to the hydrodynamic 
interactions [11], and the diffusion constant D becomes 
larger than the value in Eq. (17). This underestimation of 
D is seen at small At in the inset of Fig. 1(a). 

Since the kinetic contribution to the kinematic vis- 
cosity ?7kin/(0 is roughly proportional to D, the relation 
"Hkin "C rjcoi + ?7 P ot at large friction coefficient 7* or poten- 
tial strength a yields large Schmidt numbers Sc in DPD. 
On the other hand, small 7* and a gives Sc < 1, e.g. 
Sc = rjkin/pD = 1/2 for At <C 1 and 7/kin > ?7coi + 7? P ot- 
Sufficiently large Sc yields hydrodynamic behavior. For 
example, a large Schmidt number is required in polymer 
simulations to produce Zimm dynamics [29] — where the 
relaxation time t p of a mode with mode number p is ex- 
pected to scale (N m /p) 3 / 2 — with moderate chain 
lengths N m for an ideal chain, as demonstrated in MPC 
simulations with Sc ~ 10 [20]. Zimm dynamics was also 
reported from DPD simulations with the most typical pa- 
rameters n* = 3, a = 25 and 7* = 4.5 [8] or 5.6 [7]; 
however, the variation of the Zimm exponent with temper- 
ature observed in Ref. [8] seems to indicate that the sim- 
ulations were performed in the region between the Rouse 
and Zimm regimes. From our results above, we obtain 
Sc = 1.7 and 77 co i + 77 pot ~ 2rj kin at n* = 3, a = 25, 
7* = 4.5 with St = 0.01 and At = 0.1. Thus, this param- 
eter set is indeed in the crossover region between gas-like 
and liquid-like behavior. 

Other Thermostats. — To simulate the hydrody- 
namic behavior of complex fluids, dimensionless hydro- 
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Fig. 4: (Color online) Dependence of the diffusion constant 
D on the strength a of the soft potential (3), for n* = 3, 
St* = 0.01, and At* = 0.1. The velocity rescaling is performed 
with time step At* and cell size l c = r cut . In comparison, the 
diffusion constant D of DPD with 7* = 4.5 is also shown. 

dynamic quantities, such as the Reynolds number and 
the Schmidt number Sc, typically have to be adjusted to 
match experimental conditions. To study low- Reynolds- 
number flows of soft matter and biological systems, high 
viscosity is often required. On the other hand, DPD simu- 
lations are also often employed to study equilibrium prop- 
erties. In this case, faster diffusion and lower viscosity is 
advantageous, since it provides faster relaxation into the 
equilibrium state. Recently, a Nose-Hoover-type thermo- 
stat for the relative velocities of neighbor pairs was pro- 
posed [30,31], where the momentum is locally conserved. 
Its main idea is to thermostat systems, but to less disturb 
the original hydrodynamic transport properties (in the 
absence of any thermostat). However, the Nose-Hoover 
thermostat usually has to be combined with another ther- 
mostat to keep the temperature constant, when a sys- 
tem includes a potential with strong C 2 discontinuity like 
Eq. (3). 

The scaling of velocities [25] is an easy way to control the 
temperature in MD simulations in thermal equilibrium. In 
order to retain hydrodynamic properties, e.g. under flow, 
the main issue is momentum conservation, i. e. how to de- 
termine the velocity of the local rest frame. We suggest to 
employ the velocity scaling of the MPC method [17], which 
can be used independent of the MPC collision procedure. 
The particles are sorted into the cells of a cubic lattice with 
lattice constant l c , and the local flow velocity is identified 
with the velocity of the center of mass of all parti- 
cles in a cell. Then, the relative velocities Ui = v, — vj? 
are rescaled as Ui — * Ui \J d(N — N c )kBT/mJ2i u i 2 i where 
N is the total number of particles and iV c is the num- 
ber of cells occupied by particles. The cells are randomly 
shifted before each scaling step to ensure Galilean invari- 
ance [18]. The velocity scaling gives faster diffusion than 
DPD as shown in Fig. 4. This is particularly important 
for solvents with Lennard- Jones-type interactions, where 
the frictional contributions of a DPD thermostat adds up 
with an already high viscosity in classical MD. Velocity 
rescaling can produce temperature gradients in flow due 
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to a locally inhomogcneous energy dissipation. To reduce 
these gradients, a local version of rescaling procedure can 
be employed. To do so, many cells are grouped into larger 
bins, typically arranged sequentially in layers or columns, 
and the rescaling is performed individually for each bin. 
Alternatively, gradients on the cell scale can be avoided 
by velocity scaling with a Monte Carlo scheme [21], where 
the scaling factor fluctuates stochastically in each cell to 
reproduce the correct kinetic energy distributions. 

Summary. — We have studied the viscosity of DPD 
with finite time step, both analytically and numerically. 
The analytical results agree very well with the simulation 
data. Our theoretical results for the viscosity can be gen- 
eralized straightforwardly to other DPD methods, such as 
DPD with a multibody thermostat [22]. Thus, we have 
shown that by varying the time step At and the friction 
coefficient 7, the dynamic properties of a DPD solvent can 
be tuned, while thermodynamic properties remain unaf- 
fected. 

Furthermore, we have shown that the velocity rescal- 
ing method, which is routinely employed in MPC, can be 
adapted to MD simulations. It respects Galilean invari- 
ance and disturbs the original hydrodynamics much less 
than a DPD thermostat. 

* * * 
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